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We perform simulations of general relativistic rotating stellar core collapse and compute the 
gravitational waves (GWs) emitted in the core bounce phase of three representative models via 
multiple techniques. The simplest technique, the quadrupole formula (QF), estimates the GW 
content in the spacetime from the mass quadrupole tensor only. It is strictly valid only in the weak- 
field and slow-motion approximation. For the first time, we apply GW extraction methods in core 
collapse that are fully curvature-based and valid for strongly radiating and highly relativistic sources. 
These techniques are not restricted to weak-field and slow-motion assumptions. We employ three 
extraction methods computing (i) the Newman-Penrose (NP) scalar ^4, (ii) Regge-Wheeler-Zerilli- 
Moncrief (RWZM) master functions, and (iii) Cauchy-Characteristic Extraction (CCE) allowing for 
the extraction of GWs at future null infinity, where the spacetime is asymptotically fiat and the 
GW content is unambiguously defined. The latter technique is the only one not suffering from 
residual gauge and finite-radius effects. All curvature-based methods suffer from strong non-linear 
drifts. We employ the fixed-frequency integration technique as a high-pass waveform filter. Using 
the CCE results as a benchmark, we find that finite-radius NP extraction yields results that agree 
nearly perfectly in phase, but differ in amplitude by ~ 1 — 7% at core bounce, depending on the 
model. RWZM waveforms, while in general agreeing in phase, contain spurious high-frequency noise 
of comparable amplitudes to those of the relatively weak GWs emitted in core collapse. We also find 
remarkably good agreement of the waveforms obtained from the QF with those obtained from CCE. 
The results from QF agree very well in phase and systematically underpredict peak amplitudes by 
~ 5—11%, which is comparable to the NP results and is certainly within the uncertainties associated 
with core collapse physics. 

PACS numbers: 04.25.D-, 04.30. Db, 97.60. Bw, 02.70.Bf, 02.70.Hm 



I. INTRODUCTION 

Massive stars (Af > 8 — 10 Mq) end their nuclear burn- 
ing lives with a core composed primarily of iron-group nu- 
clei embedded in an onion-skin structure of progressively 
lighter elements. Energy generation has ceased in such 
a star's high-density core and relativistically-degenerate 
electrons provide pressure support against gravity. Sili- 
con shell burning, neutrino cooling, and deleptonization 
eventually push the core over its effective Chandrasekhar 
mass. Radial instability sets in, leading to core collapse, 
accelerated by electron capture and photodisintegration 
of iron-group nuclei (see, e.g., 

The collapsing iron core separates into a subsonically 
collapsing homologous {v oc r) inner core and superson- 
ically infalling outer core. When the former reaches nu- 
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clear density, the nuclear equation of state (EOS) stiff- 
ens, dramatically increasing central pressure support and 
stabilizing the inner core, which, due to its large inertia, 
overshoots its new equilibrium and then rebounds into 
the still collapsing outer core, launching the hydrody- 
namic supernova shock. The acceleration experienced by 
the inner core in this core bounce is tremendous, leading 
to the reversal of the collapse velocities of order 0.1c of 
its ^ 0.5 M0 of material on a millisecond timescale. 

It was realized early on that the large accelerations 
encountered in stellar collapse in combination with a 
source of quadrupole (or higher) order asphericity lead 
to the emission of a burst of gravitational waves (GWs; 
see [3: for a historical overview). Rotation, centrifugally 
deforming the inner core to oblate shape, is an obvi- 
ous source of such quadrupole asymmetry and rotating 
core collapse and bounce is the most extensively studied 
GW emission process in stellar collapse (see, e.g., [1H5] 
for recent studies and references therein). Alternatively, 
asymmetries in collapse may arise from perturbations, 
e.g., due to large convective plumes in the final phase 
of core nuclear burning, and may lead to GW emission 
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at bounce and/or seed GW-emitting prompt postbounce 
convection [3l [TOl [11] . A multitude of GW emission pro- 
cesses may be active in the postbounce, pre-explosion 
phase. These include convection/turbulence in the pro- 
toneutron star and in the postshock region, nonaxisym- 
metric rotational instabilities of the protoneutron star, 
protoneutron star pulsations, instabilities of the stand- 
ing accretion shock, and asymmetric emission of neutri- 
nos fpi and references therein). 

Of the entire ensemble of potential GW emission pro- 
cesses in stellar collapse, rotating core collapse and 
bounce is arguably the simplest and yields the clean- 
est signal, depending only on rotation, on the nuclear 
EOS, and on the mass of the inner core at bounce [5]. 
Moreover, 3D studies have shown that collapsing iron 
cores with rotation rates in the range of what is physi- 
cally plausible stay axisymmetric throughout the collapse 
phase and develop nonaxisymmetric dynamics only after 
bounce [H E] • Hence, the GW signal of rotating core 
collapse and bounce is linearly polarized and axisymmet- 
ric (2D) simulations are sufficient for its prediction. Un- 
like postbounce dynamics involving large scale and small 
scale fluid instabilities of stochastic nature, the GW sig- 
nal of rotating collapse and bounce can, in principle, be 
predicted exactly for a given set of initial data. Hence, 
it has the potential of being used in GW searches using 
matched-filtering techniques (e.g., [T^) or alternative ap- 
proaches also taking into account detailed signal predic- 
tions [UlEO]. 

Much progress has been made in recent years in the 
modeling of rotating core collapse and its GW signa- 
ture. State-of-the-art simulations are general relativistic 
(GR) [31 [3 [III [2lM7] and some studies include magnetic 
fields [Ml [Ml HE] or finite-temperature EOS, deleptoniza- 
tion, and progenitors from stellar evolutionary calcula- 
tions 21 \hj -■"> 27, ■ These improvements in the physics 
included in core collapse models provide for a more ac- 
curate and reliable dynamics underlying the emission of 
GWs. The calculation of the GW signal itself, however, is 
still being carried out predominantly in the slow-motion, 
weak- field quadrupole approximation (e.g., [5S]) that is 
of questionable quality, given the extreme densities and 
velocities involved in core collapse. The quadrupole for- 
mula (QF) "extracts" GWs based on matter dynamics 
alone, is not invariant under general relativistic gauge 
transformations, treats the emission region as a point 
source, and suffers from the fact that the definition of 
the generalized mass quadrupole moment is not unique 
in GR. 

In GR, the GW content of a spacetime can be ex- 
tracted by means of the perturbative Regge-Wheeler- 
Zerilli-Moncrief (RWZM) formalism [30H33] which is 
gauge invariant to first order or via the Newman-Penrose 
(NP) scalars approach [331 [33] which depends on the non- 
unique choice of the tetrad in which the Newman-Penrose 
scalars are evaluated. For reliable results, both RWZM 
and NP require extraction in the wave zone p9j at co- 
ordinate radii many wavelengths from the source, but 



even there, coordinate ambiguities exist. The latter are 
removed only when GWs are extracted at future null in- 
finity {J^, see [331 [3S]), where space is asymptotically 
fiat. 

Shibata & Sekiguchi |36j have used simulations of an 
oscillating polytropic neutron star model to compare QF 
and finite-radius RWZM results. For the same basic 
system, Baiotti et al. ^37] compared QF, finite-radius 
RWZM, and finite-radius NP GW extraction with each 
other and with results from a ID perturbation analy- 
sis. Both studies found that in the context of neutron 
star oscillations, the phase of the waveforms obtained 
with the quadrupole approximation agrees exceptionally 
well with that of the RWZM and NP extraction meth- 
ods. Shibata & Sekiguchi, using their particular choice of 
the generalized quadrupole moment, found a systematic 
^ 20% underprediction of the GW amplitudes by the QF. 
Baiotti ct al. |37| . who studied multiple incarnations of 
the QF, found either underprediction or overprediction of 
the amplitude, both by up to ^ 60%, depending on the 
particular choice of QF. Nagar et al. ^38J studied the per- 
formance of RWZM and QF-based GW extraction from 
oscillating polytropic tori and found qualitatively sim- 
ilar results, and quantitative differences in amplitudes 
and integrated emitted energies -Eqw between ^ 2% and 
^ 25%, again depending on the choice of quadrupole mo- 
ment. 

RWZM and NP GW extraction and comparisons with 
the QF approximation for GWs emitted in core col- 
lapse spacetimes have proven difficult. On the one hand, 
the emitted GWs are weak: Typical strain amplitudes 
are Dh ^ 10 — 1000 cm, where D is the distance to 
the source, and typical emitted energies are of order 
IQ-^o _ iQ-8 MqC^ [3], many orders of magnitude lower 
than what is expected, for example, from double neu- 
tron star coalescence |39j or binary black hole mergers 
[40l l4l] . On the other hand, the GWs have typical 
frequencies of 100 — 1000 Hz and corresponding wave- 
lengths of 300 — 3000 km, hence require extraction at 
large coordinate radii where the grid resolution of core 
collapse simulations is typically too low to allow extrac- 
tion of the relatively low-amplitude GWs emitted in core 
collapse (see, e.g., the discussion in [32]) • Shibata & 
Sekiguchi, in ^7\, were able to extract GWs with the 
RWZM formalism from an extreme core collapse model 
that developed a rotationally-induced large-scale non- 
axisymmetric deformation after bounce, emitting GWs 
with Dh ^ 20000 cm. For this model, they found that 
the QF accurately predicts the GW phase, but under- 
estimates the strain amplitude by ~ 10%. Due to the 
aforementioned difficulties, these authors were unable to 
compare RWZM with QF for more moderate, axisymmet- 
ric models. Cerda-Duran et al. [33] performed core col- 
lapse simulations using a second-order post-Newtonian 
(2PN) extension of the conformal-fiatness approximation 
to GR. Exploiting an approximate relationship of the 
non-conformal 2PN part of the metric to its GW part 
[33] , they were able to extract GWs from their 2PN met- 
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ric in standard axisymmetric rotating core collapse mod- 
els. They found very close agreement (to a few percent in 
strain amplitude) between QF and 2PN GW signals for 
almost all considered collapse models. Siebel et al. [H] 
performed nonrotating axisymmetric core collapse simu- 
lations by employing evolutions based on a fully general 
relativistic null cone formalism. They added nonspher- 
ical perturbations to the star, leading to the emission 
of GWs which they were able to extract with the Bondi 
news function at J'~^ . Comparisons to the QF suggested a 
significant discrepancy in amplitude and frequency from 
the more reliable Bondi news result. 

The results of Shibata & Sekiguchi [17] and of Cerda- 
Duran et al. |43] provide some handle on the perfor- 
mance of the QF approximation in core collapse space- 
times. The former study, while being performed in full 
GR, considered only a single extreme model. In addition, 
the authors were forced to extract GWs with RWZM at 
too small radii for completely reliable results. The latter 
study, while considering a broader ensemble of models, 
was restricted to 2PN without considering full GR, leav- 
ing room for doubts about the quality of their GW ex- 
traction technique. Finally, the results of Siebel et al. [U] 
were limited to axisymmetry without rotation and are 
unreliable in the presence of strong shocks [H] , 

In this study, we readdress GW extraction from ro- 
tating core collapse spacetimes. We perform 3-1-1 GR 
hydrodynamics simulations of rotating core collapse, for 
the first time in the core collapse context extracting GWs 
with RWZM, NP, and multiple QFs and comparing the 
results of these methods. In addition, and also for the 
first time in the present context, we utilize the Cauchy- 
Gharacteristic Extraction (CCE) approach that 
propagates the GW information to J'^ for completely 
gauge independent and unambiguous GW extraction. 

In choosing our models set, we are guided by Cerda- 
Duran et al. [43], and draw precollapse configurations 
from the set of [21]. These models are GR n = 3- 
polytropic iron cores in rotational equilibrium and we 
evolve them with an analytic hybrid polytropic/F-law 
EOS used in many previous studies of rotating core col- 
lapse [13 [m m m 113 [SD] . For physically accurate 
GW signal predictions to be used in GW data analy- 
sis, a microphysically more complete treatment is war- 
ranted. Fortunately, recent results of studies employing 
such modehng technology (e.g., [H [51 1251 ^Jl 1ST]) show 
that, with a proper choice of EOS parameters, hybrid- 
EOS models are able to qualitatively and to some ex- 
tent quantitatively reproduce the GW signals obtained 
with the much more complex and computationally inten- 
sive microphysical studies. Hence, for the purpose of this 
study, we resort to the simpler hybrid-EOS models. 

Our simulations employ the open-source ZelmEini GR 
core collapse simulation package [35] that is based on 
the Cactus Computational Toolkit |53l I54j and the 
Einstein Toolkit [55]. While using the fuh 3-1-1 GR 
formalism, we limit our simulations to an octant of the 
3D cube, using periodic boundary conditions on two of 



the inner faces of the octant and reflective boundary con- 
ditions on the third face. This limits 3D structure to even 
£ and m that are multiples of 4, which is not a limitation 
for the current study, since rotating core collapse and 
the very early postbounce evolution are likely to proceed 
nearly axisymmetrically [4l[6l[56]. We note that, even 
though the GW signal in rotating core collapse is domi- 
nated by the {£ = 2,m = 0) polarization mode, there 
is no reason to expect different behavior for other GW 
multipoles or polarizations and our results should trans- 
late to the non-axisymmetric case. 

The results of our simulations indicate that NP extrac- 
tion yields results that agree well with those obtained 
from the most sophisticated GGE method. We observe 
differences in amplitude of 1 — 7%, depending on the 
model, while the agreement in phase is nearly perfect. 
We also find that the RWZM formalism yields unphysical 
high-frequency signal components that make this method 
less suitable for core collapse simulations where the sig- 
nal is very weak. Finally, we note that the quadrupole 
approximation yields surprisingly close results to those 
obtained from CCE. While the phases nearly perfectly 
agree, the amplitude shows differences of 5 — 11%. 

This paper is structured as follows. In Sec. |ll| we dis- 
cuss our methodology, initial data, and EOS details. Sec- 
tion mil discusses the various GW extraction methods 
that we employ. In Sec. |IV[ we present our results and 
discuss them in detail. Finally, in Sec. |V| we summarize 
and review our findings. 



II. METHODS 

We adopt the ADM 3+1 fohation of spacetime [57]. AU 
equations assume c — G = 1 unless noted otherwise. In 
the following, Latin indices run from 1 to 3 while Greek 
ones run from to 3. We adhere to abstract index nota- 
tion, g^^ is the 4-metric, jij is the 3-metric and Kij the 
extrinsic curvature. 



A. Infrastructure and Mesh Refinement: Cactus 
and Carpet 

Our code uses the Cactus Computational Toolkit 
[531 I54j to manage the complexity inherent in large 
software projects. Cactus is an open source high- 
performance computing environment designed for scien- 
tists and engineers; its modular structure enables paral- 
lel computation across different architectures, and facil- 
itates collaborative code development between different 
groups. Indeed, our code uses a set of components of the 
public EinsteinToolkit ^55. ,58J, a community project 
developing and supporting open software for relativistic 
astrophysics, such as e.g. the curvature and hydrodynam- 
ics evolution methods described below. Many improve- 
ments made in the course of the research for this paper 
were contributed back to the community. 
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In particular, Cactus allows us to clearly separate be- 
tween physics components and computational compo- 
nents in our code. Distributed memory parallelism in 
Cactus is provided by a driver component which im- 
plements the data structures discretising the manifold 
on which the computational state vector lives. In our 
case, this is the Carpet driver [59H61) providing block- 
structured adaptive mesh refinement (AMR) and multi- 
block discretization. Carpet parallelizes using a hybrid 
approach combining MPI and OpenMP, where inter-node 
communication is handled via MPI and intra-node com- 
munication via OpenMP or via MPI, depending on the 
particular system and on details of the simulation setup. 

Carpet implements Berger-Oliger style AMR [62] . 
where the fine grids are aligned with coarse grids, refined 
by factors of two. Carpet also implements subcycling 
in time, where finer grids take two time steps for every 
coarse grid step. The latter greatly improves efficiency, 
but also introduces significant complexity into the time 
evolution method. The refined regions can be chosen 
and modified arbitrarily, which we use here to add addi- 
tional, finer levels during evolution as successively higher 
resolutions are required to capture the collapse. This is 
described in more detail in |42j . 

We use fifth-order accurate spatial interpolation for 
spacetime variables and third order essentially non- 
oscillatory interpolation for hydrodynamics variables. 
Time interpolation, which is necessary to provide bound- 
ary conditions to fine levels at times where there is no 
coarse level, is second-order accurate. We apply no time 
refinement between levels 3 and 4, which corresponds to 



reducing the Courant-Friedrichs-Lewy factor on levels 3 
and coarser by a factor of 2. This increases the accuracy 
on level 3 where we extract gravitational waves. In to- 
tal, we use 9 refinement levels (including the base grid), 
an outer boundary radius of 384OA'/0 (^ 5700 km) and 
a finest zone size of O.25M0 (^ 370 m) in our baseline 
resolution. 



B. Curvature Evolution: McLachlan 

1. Evolution System 

We evolve the full Einstein equations in a 3 -|- 1 split 
(a Cauchy initial boundary value problem) [33], using 
the BSSN formulation [64], a 1 log shcing [65], and F- 
drivcr shift condition [65] . This leads to the following set 
of evolved variables: 



log 



det7ij 



K 

Ai-j :— 



(1) 

(2) 
(3) 

(4) 

, (5) 

Our exact evolution equations are as described by 
Eqs. (3) to (10) of which we list here for complete- 



K 



doK 

do(j) 
do A J 



D'D^a + 2^^(f> ■ D'a 



+ a A'^A, 



6 6 

~2aAj + 27fe(,a,)/3'= - '^^rjdkP'' 



aRij + CiR'lj — DiDjU + 49(^0 • Dj^a 

2 r 



TF 



2A^dja + 2a 



(m - l)dkA'' 



ki 



9m 

— D'i^ + m{T\iA''-^ + GA'^d, 



(6) 
(7) 

(8) 
(9) 

(10) 
(11) 



(12) 



(13) 



r 



with the momentum constraint damping constant set to 
771 = 1. The stress energy tensor T^i^ is incorporated via 



the projections 





- 2P'To^ + P'P^T') 


(14) 


S := g^'T^j 




(15) 


--(To, 


- P'T,,) . 


(16) 
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We have introduced the notation do = dt — ji-'dj. All 
quantities with a tilde ~ refer to the conformal 3-metric 
7ij, which is used to raise and lower indices. In par- 
ticular, Di and F*^ refer to the covariant derivative and 



the Christoffel symbols with respect to 7ij. The expres- 
sion [• • • J"'"^ denotes the trace- free part of the expression 
inside the parentheses, and we define the Ricci tensor 
contributions 



kij 



(17) 
(18) 



This is a so-called (/)- variant of BSSN. The evolved gauge 
variables are lapse a, shift /3*, and a quantity related 
to the time derivative of the shift. The gauge parameters 
/, G, H, and rj are determined by our choice of a 1 + log 
shcing: 



f{a,<l),xn := 2/a 
Koixf") 

and F-driver shift condition: 



ri{B\a,x'') 



(3/4) 
exp{4(/)} 



(19) 
(20) 



(21) 
(22) 
(23) 



The expression q{r) attenuates the F-driver depending 
on the radius as described below. 

The F-driver shift condition is symmetry-seeking, driv- 
ing the shift /3* to a state that renders the conformal con- 
nection functions F* stationary. Of course, such a station- 
ary state cannot be achieved while the metric is evolving, 
but in a stationary spacetime the time evolution of the 
shift and thus that of the spatial coordinates a;* will be 
exponentially damped. This damping time scale is set by 
the gauge parameter 77 (see Eq. 23 1 which has dimension 
1/T (inverse time). As described, e.g., in [67l[68], this 
time scale may need to be adapted in different regions of 
the domain to avoid spurious high-frequency behavior in 
regions that otherwise evolve only very slowly, e.g., far 
away from the source. 

Here we use the simple damping mechanism described 
in Eq. (12) of 68 , which is defined as 



q{r) 



j 1 for r < R (near the origin) 
) R/r for r > R (far away) 



(24) 



with a constant R defining the transition radius between 
the interior, where g ~ 1, and the exterior, where q falls 
off as 1/r. Eq. [23] describes how q appears in the gauge 
parameters. In this paper we use R = 250 Af© (i? — 
369.2 km). 

We implement the above BSSN equations and gauge 
conditions in the McLachlan code [BHl HI] which is freely 
available as part of the EinsteinToolkit. McLachlan is 



auto-generated from the definition of the variables and 
equations in the Mathematica format by the Kranc code 
generator [70H72] . Krcoic is a suite of Mathematica pack- 
ages comprising a computer algebra toolbox for numeri- 
cal relativists. Kranc can be used as a "rapid prototyp- 
ing" system for physicists or mathematicians handling 
complex systems of partial differential equations, and 
through integration into the Cactus framework one can 
also produce efficient production codes. 

We use fourth-order accurate finite differencing for the 
spacetime variables and add a fifth-order Kreiss-Oliger 
dissipation term to remove high frequency noise. We 
use a fourth-order Runge-Kutta time integrator for all 
evolved variables. 



2. Initial Conditions 

We set up our initial condition from the ADM variables 
gij, Kij, lapse a, and shift /?*, as provided by the initial 
data discussed in Sec. IIIDl From these we calculate the 
BSSN quantities via their definition, setting = 0, and 
using cubic extrapolation for F* at the outer boundary. 
This extrapolation is necessary since the F' are calcu- 
lated from derivatives of the metric, and one cannot use 
centered finite differencing stencils near the outer bound- 
ary. We assume that one could instead also use one-sided 
derivatives to calculate F* on the boundary. 

The extrapolation stencils distinguish between points 
on the faces, edges, and corners of the grid. Points on the 
faces are extrapolated via stencils perpendicular to that 
face, while points on the edges and corners are extrapo- 
lated with stencils aligned with the average of the nor- 
mals of the adjoining faces. For example, points on the 
{+x, +y) edge are extrapolated in the (1, 1,0) direction, 
while points in the (-fx, +y + z) corner are extrapolated 
in the (1, 1, 1) direction. Since several layers of boundary 
points have to be filled for higher order schemes (e.g., 
three layers for a fourth order scheme), we proceed out- 
wards starting from the innermost layer. Each subse- 
quent layer is then defined via the points in the interior 
and the previously calculated layers. 
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3. Boundary Conditions 

During time evolution, we apply a Sommerfeld-type 
radiative boundary condition to all components of the 
evolved BSSN variables as described in The main 

feature of this boundary condition is that it assumes ap- 
proximate spherical symmetry of the solution, while ap- 
plying the actual boundary condition on the boundary 
of a cubic grid where the face normals are not aligned 
with the radial direction. This boundary condition de- 
fines the right hand side of the BSSN state vector on the 
outer boundary, which is then integrated in time as well, 
so that the boundary and interior are calculated with the 
same order of accuracy. 

The main part of the boundary condition assumes that 
we have an outgoing radial wave with some speed Vq: 



X 



u(r - vot) 



(25) 



where X is any of the tensor components of evolved vari- 
ables, Xq the value at infinity, and u a spherically sym- 
metric perturbation. Both Xq and vq depend on the par- 
ticular variable and have to be specified. This implies 
the following differential equation: 



dtX 



-v'd.X - Vq 



X-Xn 



(26) 



where — VQX^/r. The spatial derivatives di are eval- 
uated using centered finite differencing where possible, 
and one-sided finite differencing elsewhere. We use sec- 
ond order stencils in our implementation. 

In addition to this main part, we also account for those 
parts of the solution that do not behave as a pure wave, 
e.g.. Coulomb type terms caused by infall of the coor- 
dinate lines. We assume that these parts decay with a 
certain power p of the radius. We implement this by con- 
sidering the radial derivative of the source term above, 
and extrapolating according to this power-law decay. 

Given a source term (dtX), we define the corrected 
source term (dtX)* via 



(dtxy = {dtX) + 



r — n^dir 



n'd,{dtX), (27) 



where rf is the normal vector of the corresponding 
boundary face. The spatial derivatives di are evaluated 
by comparing neighbouring grid points, corresponding to 
a second-order stencil evaluated in the middle between 
the two neighbouring grid points. We assume a second- 
order decay, i.e., we choose p = 2. 

As with the initial conditions above, this boundary 
condition is evaluated on several layers of grid points, 
starting from the innermost layer. Both the extrapo- 
lation and radiative boundary condition algorithms are 
implemented in the publicly available NewRad component 
of the Einstein Toolkit. 

This boundary condition is only a coarse approxima- 
tion of the actual decay behavior of the BSSN state vec- 
tor, and it does not capture the correct behavior of the 



evolved variables. However, we observe that this bound- 
ary condition leads to stable evolutions if applied suf- 
ficiently far from the source. Errors introduced at the 
boundary (both errors in the geometry and constraint vi- 
olations) propagate inwards with the speed of light |66j . 
Gauge changes introduced by the boundary condition, 
which are physically not observable, propagate faster, 
with a speed up to for our gauge conditions. 



C. General-Relativistic Hydrodynamics: GRHydro 

We employ the open-source GR hydrodynamics code 
GRHydro that is part of the EinsteinToolkit 55J and is 
an updated version of the code Whisky described in |73j . 

The equations of ideal GR hydrodynamics evolved by 
GRHydro are derived from the local GR conservation laws 
of mass and energy-momentum. 



J'' = 0, 



V^Tf"" = 



(28) 



where Vu denotes the covariant derivative with respect 



to the 4-metric. 
the 4- velocity 
phu^ 



pu 



^ is the mass current with 



and the rest-mass density p. T^'^ 



{'Pg^'^ is the stress-energy tensor. The quantity 
h — 1 4- e -f P/ p is the specific enthalpy, P is the fluid 
pressure and e is the specific internal energy. 

We choose a definition of the 3-velocity that corre- 
sponds to the velocity seen by an Eulerian observer at 
rest in the current spatial 3-hypersurface |7lj , 



(29) 



where W ^ {I — v^Vi)^^^^ is the Lorentz factor. In 
terms of the 3-velocity, the contravariant 4-velocity is 
then given by 



W 

a 



u' = W [v 



a 



and the covariant 4-velocity is 

uq — Wiv'-Pi — a) , Ui = Wvi 



(30) 



(31) 



The GRHydro scheme is written in a first-order hyper- 
bolic flux-conservative evolution system for the conserved 
variables D, 5*, and f in terms of the primitive variables 



D = ^pW, 
S' = ^phW^v\ 
T = ^{phW^-P)-D, 



(32) 



where 7 is the determinant of 7ij . The evolution system 
then becomes 



(33) 
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with 



U 

pi 



dxf^ 



(34) 



Here, = — /3''/a and are the 4-Christoffel sym- 
bols. The above equations are solved in semi-discrete 
fashion. The spatial discretization is performed by means 
of a high-resolution shock-capturing (HRSC) scheme em- 
ploying a second-order accurate finite- volume discretiza- 
tion. We make use of the Marquina flux formula for the 
local Riemann problems and piecewise-parabolic cell in- 
terface reconstruction (PPM) . For a review of such meth- 
ods in the GR context, see [7S]. The time integration and 
coupling with curvature are carried out with the Method 
of Lines [7S] . 



D. Equation of State and Initial Stellar Models 

For the purpose of this study, we employ the simple an- 
alytic hybrid EOS [ST] that combines a 2-piece piece- 
wise polytropic pressure Pp with a thermal component 
Pth, i.e., P — Pp + Pth- To model the stiffening of the 
EOS at nuclear density p^uc = 2x 10^^ gcm~'^, we assume 
that the polytropic index 7 jumps from 71 below nuclear 
density to 72 above. As detailed in [TH], it is possible to 
construct an EOS that is continuous at Pnuc, 



P 



1 - 7th 



7- _ 
+ (7th 



Kp. 



nuc r 



(7th - l)(7-7i) 
(7i-l)(72-l) 



Kp 



nuc P 



l)pe. 



(35) 



In this, e = ep -|- eth denotes the total specific internal 
energy which consists of a polytropic and a thermal con- 
tribution. K = 4.897 X 10"'^^ [cgs] is the polytropic con- 
stant for a polytrope of relativistic degenerate electrons 
at Ye = 0.5 . The thermal index 7th = 1.5 corresponds 
to a mixture of relativistic (7 = 4/3) and non-relativistic 
(7 = 5/3) gas. This EOS mimics the effects of the stiff- 
ening of the physical EOS at p-^-ac and can handle the 
significant thermal pressure contribution introduced by 
shock heating in the postbounce phase. Provided appro- 
priate choices of EOS parameters (e.g., [17]), the hybrid 
EOS leads to qualitatively correct collapse and bounce 
dynamics. Consequently, this leads to GW signals that 
are similar in morphology, characteristic frequencies and 
amplitudes to those computed from more microphysically 
complete simulations [H [51 [HI] . 

We employ n = 3 (7ini = 7 = 4/3) polytropes in rota- 
tional equilibrium generated via Hachisu's self-consistent 
field method [79 , 80 1 that not only provides fluid, but also 
spacetime curvature initial data. The polytropes are set 



TABLE L Initial parameters of differentially rotating stellar 
cores used for the core collapse simulations. The models are 
described by three quantities: the degree of differential ro- 
tation A, the ratio T/|M^| of rotational to potential energy, 
and the sub- nuclear adiabatic index 71 during the collapse. 
For convenience we also report the wave-signature type of the 
three models and the mass M present on the computational 
grid. 



Model 


Type 


A [10^ km] 


T/\W\[%] 


71 


M[Mq] 


A1B3G3 


I (weak) 


50.0 


0.89 


1.31 


1.46 


A1B3G5 


III 


50.0 


0.89 


1.28 


1.46 


A3B3G3 


I (strong) 


0.5 


0.89 


1.31 


1.46 



up with the rotation law discussed in [3TJ [SD] and are 
parametrized via the differential rotation parameter A 
and the initial ratio T/|W^| of rotational kinetic energy 
T to gravitational binding energy \W\. While being set 
up as marginally stable polytropes with 7ini = 4/3, dur- 
ing evolution, the initial sub-nuclear polytropic index 71 
is reduced to 71 < 7ini to accelerate collapse. Follow- 
ing previous studies [HJ [551 [iO], we use 72 — 2.5 in the 
super-nuclear regime. 

From the initial stellar configurations of [211 [S3 we 
draw a subset of three models that cover the range of 
astrophysically expected GW signals from rotating iron 
core collapse [5] and accretion- induced collapse [51 . Our 
choices have been used previously in a comparison study 
of full GR and conformally-flat simulations |25j : 

• Model A1B3G3 is in near uniform rotation with 
A = mx 10^ km, has T/\W\ = 0.9%, and, once 
mapped to the evolution grid, uses a sub-nuclear 
adiabatic index 71 = 1.31. Its GW signal is of the 
standard "Type-I" morphology [HI [23 [50] and of 
moderate strength (see [211 (25] for details) . 

• Model A3B3G3 also uses 71 = 1.31. It is strongly 
differentially rotating, with its initial central angu- 
lar velocity dropping by a factor of two over A = 
500 km. This, in combination with T/\W\ = 0.9%, 
leads to rapid rotation in the inner core, resulting 
in a very strong GW signal at core bounce and dy- 
namics that are significantly affected by centrifugal 
effects. It produces a "Type-I" GW signal with a 
centrifugally-widened broad peak at core bounce. 

• Model A1B3G5 has the same rotational setup as 
model A1B3G3, but its sub-nuclear adiabatic in- 
dex is reduced to 71 = 1.28. This leads to rapid 
collapse, to a very small inner core at core bounce, 
and to a weak "Type-Ill" GW signal [H [50] akin 
to that potentially emitted by an accretion-induced 
collapse event |5TI . 

For convenience, key model parameters are summarized 
in Table m 
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III. GRAVITATIONAL WAVE EXTRACTION 
A. The Quadrupole Approximation 

The quadrupole approximation is the only means of 
extracting GWs in Newtonian or conformally-flat GR 
simulations, but has found wide application also in GR 
simulations of stellar collapse [H [T71 |23J US] . 

The coordinate-dependent quadrupole formula esti- 
mates the GW strain seen by an asymptotic observer 
by considering exclusively the quadrupole stress-energy 
source. It neglects any non-linear GR effects. This 
approximation is valid strictly only in the weak-field 
^ <C 1 and slow-motion f <C 1 limit [29] where space- 
time is essentially flat. 

The quadrupole formula is given in the transverse- 
traceless (TT) gauge by 



Ijk{t-R/c) 



TT 



where 



Ijk / pit, x) 



(fa 



(36) 



(37) 



denotes the reduced mass-quadrupole tensor. Since we 
are working in the weak-field, slow-motion approxima- 
tion, the placement of tensor indices is arbitrary. Ijk is 
not uniquely defined in GR and the choice of the density 
variable p is not obvious. Following [H [51 [5T1 [231 US] , we 
set p — y/^Wp = D, because, (i), this is the conserved 
density variable in our code, and (ii), ^J^<Px is the nat- 
ural volume element. See [37] for other potential choices 
and their relative performance for GWs from oscillating 
polytropes. 

The reduced mass-quadrupole tensor can be computed 
directly from the computed distribution I?(t,x). In or- 
der to eliminate the effects of numerical noise when dif- 



ferentiating Eq. ( 37 1 twice in time, we make use of the 



continuity equation to obtain the first time derivative of 



Eq. (37) without numerical differentiation [22l ISTl ,82 . 



dt 



\{x^i^)5^^ 



d^3 



29 



(38) 
Note that 



where we set = as defined by Eq. 
we have switched to contravariant variables in the inte- 
grand as these are the ones present in the code. Finally, 
the remaining time derivative needed for evaluating the 
quadrupole GW strain (Eq. [36[ ) is performed numerically. 

In order to assess the sensitivity of the predicted waves 
on the particular choice of the velocity variable in 



Eq. (38), we implement two modified versions. In variant 



VS, we use Shibata et al.'s definition of the 3- velocity 
(e.g., [23]) that differs from ours by a gauge term. In 
variant PV, we follow [7S] and employ physical veloc- 
ity components (individually bound to u < c) that. 



in Cartesian coordinates, are given by {vx,Vy,Vz} ~ 
{^/jiiv^ , y/j22v'^ , y/j33V^}, assumiug that the 3-metric 
is nearly diagonal (which is the case in our gauge; see 
Sec.|ra. 



B. The Regge-Wheeler-Zerilli-Moncrief Formalism 

A particular Ansatz for analyzing gravitational radi- 
ation in terms of odd and even multipoles in the far- 
field of the source was originally developed by Regge, 
Wheeler [3D] and Zerilli [ST] [35] , respectively. Moncrief 
subsequently provided a gauge-invariant reformulation 
[83] (see [84] for a review). Assuming that, at large dis- 
tances from the source, the GW content of the spacetime 
can be viewed as a linear perturbation to a fixed back- 
ground, we can write 



5^1^ = 9n 



(39) 



where g'^^ is the fixed background metric and its 
linear perturbation. The background metric 17°^ is usu- 
ally assumed to be of Minkowski or Schwarzschild form, 
which we can write as 



ds-" 



'Ndt^ + Adr^ + r^{de^ +i 



(40) 



By splitting the spacetime into timclikc and radial and 
angular parts, it is possible to decompose the metric per- 
turbation hf^i, into odd and even multipoles, i.e., we can 
write 



(o) 



(41) 



The even and odd multipole components are defined ac- 
cording to their behavior under a parity transformation 
{0,4>) —T' {tt — 9,Tr + (f)). Odd multipoles transform as 
(—1)''+^ while even multipoles transform as (—1)^. Both 
multipole components can be expanded in terms of vector 
and tensor spherical harmonics (e.g., 

Given the Hamiltonian of the perturbed Einstein equa- 
tions in ADM form ^85J , it is then possible to derive vari- 
ational principles for the odd and even-parity perturba- 
tions ^ET to give equations of motions that are similar to 
wave equations with a scattering potential. 

The solutions to the odd and even-parity wave equa- 
tions are given by the Regge- Wheeler-Moncrief and the 
Zerilli-Moncrief master functions, respectively. The odd- 
parity Regge- Wheeler-Moncrief function reads 



l2{i+iy. 1 / 2M\ 



and the even-parity Zerilli-Moncrief function reads 
/2(£ + l)! 



em 
"il 



{£-2)1 A[r(A-2) + 6M] 



(42) 



(43) 
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where A = + 1), and where 



4r 



(44) 



with 



Jm 



1 r 

2 



, (45) 



(46) 



These master functions depend entirely on the spheri- 
cal part of the metric given by the coefficients N and A, 
and the perturbation coefficients for the individual met- 
ric perturbation components {h{'^)'^°\ {hi"')^°\ (/if')(°\ 



tl2 : 



K*^"^, and G which can be 



obtained from any numerical spacetime by projecting out 
the Schwarzschild or Minkowski background [86j . For ex- 
ample, the coefhcient can be obtained via 



H. 



(47) 



where is the radial component of the numerical met- 
ric represented in the spherical-polar coordinate basis, 
are spherical harmonics, and dfl is the surface line 
element of the extraction sphere. The coefficient A 
represents the spherical part of the background metric 
and can be obtained by projection of the numerical met- 
ric component g^ on Ygo over the extraction sphere 



A^ 



1 

An 



(48) 



Similar expressions hold for the remaining perturbation 
coefficients. 

The odd- and even-parity master functions Eq. ( 42 ) 



and Eq. (43 1 can be straight-forwardly related to the 



gravitational-wave strain and are given by 



(49) 



where {0,(f>) are the spin- weight s = —2 spherical 

harmonics. We note that this relation is strictly true 
only at an infinite distance from the source. Since our 
numerical domain is finite in size, we choose some, ideally 
large, but finite radius. In Sec. |IV C[ we check how well 
the GWs extracted with the RWZM formalism asymptote 
with increasing extraction radius. 

In the present work, our models exclusively trigger the 
even-parity master function , and is zero. In this 
case, we can simplify Eq. ( 49 ) and obtain 



1 



V2i 



(50) 



relating the strain directly to 



C. Newman-Penrose Scalars 

Another method for calculating the gravitational wave- 
forms is based on the conformal structure of asymptoti- 
cally flat spacetimes as established by Bondi, Sachs and 
Penrose [321 1571 IHB] . This method is conveniently rep- 
resented in terms of spin-weighted scalars as introduced 
by Newman and Penrose |34j . In the following we refer 
to it as NP extraction. According to the peeling theo- 
rem [S7J IHB] , a certain component of the conformal Weyl 
tensor obeys the slowest l/r fall-off from the source, and 
hence is identified as outgoing gravitational radiation: 



^'4 



*1 , *0 



+ Z:^ + ^ + ZT + Z^+0{r-'). (51) 



Here, the slowest fall-off is obeyed by the NP scalar 'I'4, 
which is defined as^ 



(52) 



where Cajs-yS is the conformal Weyl tensor associated with 
the 4-metric gap and n, ffi are part of a null-tetrad |34U35| 
H., n, m, m which satisfies ~£ ■ n = 1 = m ■ ffi while all 
other inner products vanish. In addition, this tetrad is re- 
lated to the 4-nietric via gat = lanh+l(,7ia—rnarhi,—m},rha. 
At future null infinity J''^, the topology of the spacetime 



is a time succession of spheres, S 



Hence the simplest 



choice for the null tetrad at J'~^ is such that it resembles 
the unit sphere metric. Moreover, the simplest choice for 
a coordinate system at J'"'" is given by the Bondi gauge 
|87, 88 , which makes use of an areal radius coordinate. 

In most current numerical relativity simulations, the 
radiation is computed at a finite radius where the Bondi 
coordinates are usually not imposed. Rather, we use the 
gauge as evolved by the 1 -I- log slicing and P-driver con- 



ditions discussed in Sec. II B In practice, we impose a 
simple polar-spherical coordinate system with constant 
coordinate radius — + + , which does not take 
into account the background geometry, and hence does 
not make use of an areal radius. Thus, the gravitational 
radiation as computed on these coordinate spheres is not 
measured in the correct gauge, and leads to a systematic 
error that needs to be assessed. Note that it is principally 
possible to transform to the correct gauge [90] . 

In our construction of an approximate tetrad, we follow 
common practice (e.g. [9TH93]) and use a triad of spatial 
vectors u, v, w obtained via a Gram-Schmidt orthonor- 
malization starting from 



u 



[X, y, z\ 



V ~ [xz, yz, — X 



(53) 
(54) 
(55) 



^ Our definition proceeds along the lines of Appendix C of Ref. |89l 
but for comparison with the quadrupole results, we define the 
Newman-Penrose scalar with the opposite sign of their Eq. (CI). 



10 



where x, y, z are Cartesian coordinates of the computa- 
tional grid and e*™" is the three-dimensional Levi-Civita 
symbol. The tetrad is given in terms of this triad and 
the timclike normal vector by 



It is convenient to decompose the two GW polarizations 



v2 



(56) 
(57) 
(58) 



A straightforward calculation shows that we are thus 
able to express 4*4 exclusively in terms of the "3-1-1" vari- 
ables according to 

^-4 - ^ [S,„„(u;'"u;" - i;™«") + S™„(«'"w" + 7«"i;")] 
+ ^ {E^niy^'w^ - u;™«") -I- B.^,\w^w^ -f v^v^)\ (59) 

where the electric and magnetic parts of the Weyl tensor 
are defined as EH 



(60) 
(61) 



Here the * denotes the Hodge dual and J-^a = S^'a+n^n^ 
is the projection operator. The Gauss-Codazzi equations 
(see e.g. [95]) enable us to calculate the electric and mag- 
netic parts from the "3-1-1" variables according to 



7""(i^,,if„„-if™i^,„), (62) 



B 



^, ~ 7»fce''""^mi^ry. (63) 



In a given numerical simulation, we calculate ^'4 from 
Eq. ( 59 1 on a set of coordinate spheres defined by i?ox = 



const. On each of these spheres, we use spin-weighted 
spherical harmonics -2Ye.m{0, 0) of spin weight —2 in or- 
der to decompose the resulting wave signal into multi- 
poles 

*4(i, e, </>) =5^*^ (t)-2y£m(^,0), 

^^™(t) = / -^^{1, 9, ^)^2%m{e, <jy)dn. (64) 



In all our simulations, the wave signal is dominated by 
the ^ = 2, TO = mode whose angular dependence is 
given by 



15 



(65) 



The NP scalar ^'4 is related to the gravitational wave 
strain via 



into multipoles in analogy to Eq. ( 64 1 



00 £ 

t=2 rn=-e 

(67) 

These multipoles are related to those of the NP scalar by 

(68) 



Note that the final result is not fully gauge-invariant 
and contains an unknown amount of systematic error. 
The reasons are two-fold: First, we did not choose a 
proper Bondi null tetrad on our extraction spheres, and, 
second, the extraction spheres have finite radius, thus are 
neglecting non-linear back-scattering effects of the gravi- 
tational field in the wave zone. However, since our coor- 
dinate frame asymptotically approaches the Minkowski 
spacetime, both errors can be minimized by performing 
extrapolation to future null infinity using a set of 
extraction spheres at finite radii. Unfortunately, even 
if the extrapolation is accurate, an uncertain amount of 
residual error may remain. In Sec. |IV C[ we check how 
well the extracted waves approximate their asymptotic 
shape and magnitude. 



D. Cauchy-Characteristic Extraction 

To circumvent the problem of finite-radius extraction 
and to eliminate this systematic error, we apply the tech- 
nique of CCE ^ ,47n49. i96n97| to obtain the NP scalar 
^4 as discussed in the previous section^, in this case di- 
rectly evaluated at future null infinity J'^. The CCE 
technique couples an exterior characteristic evolution of 
the full Einstein equations to the interior strong-field 3+1 
Cauchy evolution of the spacetime. 

Characteristic evolutions are based on nuU- 
hypersurface foliations of spacetime and have the 
advantage of allowing for a compactification of the 
radial coordinate component, thus allowing to include 
future null infinity J^^ on the computational grid [45j . 
Unfortunately, the characteristic formulation gives rise 
to the formation of caustics, i.e., the null rays on 
which the coordinate system is based can intersect in 
strong-field regions, leading to coordinate singularities. 
The scheme is therefore not well suited for the evolution 
of the actual GW source. Characteristic evolutions, on 
the other hand, are well adapted to the far-field region of 
spacetime and can efficiently evolve the metric fields out 



*4 



(66) 



^ Alternatively during this procedure, we also compute the Bondi 
news function |34| A/", which is related to the GW strain by only 
one integration in time. We find that the news function is less 
robust when residual matter is present at the world-tube location 
and therefore restrict our attention to \E'4 only. 
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to where it is possible to obtain ^'4 (and, hence, K) 
in a mathematicaUy unambiguous and gauge-invariant 
way gainiHMl- 

We therefore proceed as foUows: we evolve the interior 
region containing the collapsing matter with the standard 
Cauchy formulation as described in Sec. |IIB| and |IIC| 
During the Cauchy evolution, we store the 3-metric com- 
ponents including lapse and shift on coordinate spheres 
with fixed radius i?r defining the world-tube F. 

This world-tube forms the inner boundary for the sub- 
sequent characteristic evolution of the Einstein equa- 
tions. The full 4-metric can be reconstructed from the 
stored 3-metric components, the lapse and the shift at 
the inner boundary. Upon construction of proper ini- 
tial data on an initial null hypersurface, which here we 
simply assume to be conformally flat, we have then fully 
specified any data necessary to evolve the fields out to 
. More details on the exact mathematical procedure 
can be found in [351 [95] . The characteristic field equa- 
tions are solved numerically using the PITT null evolution 
code [46' . The numerical implementation of CCE includ- 
ing results from binary black hole mergers is discussed in 
[37H3HI • For the characteristic computational grid, we use 
= 397 points in the radial direction and iVang — 73 
points in each angular direction for the two stereographic 
patches covering the sphere. The characteristic timestep 
size equals that of the Cauchy evolution. 

After each characteristic timestep, the NP scalar 
is evaluated directly from the metric at and trans- 
formed to the desired Bondi gauge [98]. Thus, the CCE 
method is free of gauge and near-zone effects and repre- 
sents the most rigorous extraction technique. However, 
there is still some remaining systematic error that is due 
to the presence of matter at the world-tube locations. 
Since the current set of characteristic equations does not 
take into account any form of matter contribution, a non- 
zero stress-energy tensor introduces an unknown error. 
We therefore have to perform checks of the dependence 
of the waveforms on the world-tube locations. In princi- 
ple, it is possible to also incorporate matter on the char- 
acteristic side [lOOj . which we leave to future work. 

We note that CCE does not remove the artifical outer 
grid boundary from the Cauchy evolution. Thus, incon- 
sistencies arising from this boundary can, in principle, 
still influence the interior domain. It is possible to cir- 
cumvent this problem by enlarging the computational 
domain so that the outer boundary is causally discon- 
nected from the world-tube locations (see, e.g., [5T[ in 
the context of binary black holes) . In simulations of core 
collapse, however, this is currently not computationally 
feasible, but experiments with varied outer boundary lo- 
cations have shown that boundary effects are negligible 
for our current choice of boundary location. 

Finally, we point out that inconsistencies in the char- 
acteristic and Cauchy initial data may lead to a loss of 
some non-linear effects. Even though we expect these 
problems to be very small (see [ISj)- These and the outer 
boundary issues highlighted in the above can be fully 



accounted for only by employing Cauchy-characteristic 
matching (CCM) (e.g. [45]). This technique uses the 
characteristic evolution as a generator for Cauchy bound- 
ary data at the world-tube, i.e. the world-tube becomes 
a two-way boundary between Cauchy and characteristic 
evolution. In practice, CCM has not yet been successfully 
implemented. 



E. Remarks on Integration and Physical Units 

The NP scalar 5*4™ must be integrated twice in time to 
yield the strain /i, which introduces an artificial "mem- 
ory" [lOlj . i.e., a non- linear drift of the signal so that the 
wavetrain deviates from an oscillation about zero. This 
behavior cannot be explained by the two unknown inte- 
gration constants resulting, at most, in a linear drift. 

As suggested in [102) . this non- linear drift is a conse- 



quence of random-walk-like behavior induced by numer- 
ical noise. In the present work, we make use of methods 
that are strictly valid only in pure vacuum and at our 
extraction spheres the average matter densities are non- 
zero (see Table III). This systematic error can lead to an 



additional artificial low- frequency drift. In order to elim- 
inate this effect, we use fixed-frequency integration (FFI) 
as proposed in 102 . The NP variable ^'4°(t) is Fourier 
transformed, the resulting spectrum is divided by /q for 
frequencies f < fo and divided by otherwise. An 
inverse Fourier transform then yields the strain him es- 
sentially free of spurious drifts and oscillations, given a 
proper choice of /q. 

Finally, we need to address the question of units. The 
gravitational wave strain — ih-^ is by construction di- 
mensionless. For comparison of waveforms at different 
extraction radii, it is convenient to compensate for the 
1/D fall-off of the strain, where D is the distance from 



{Im) 



and 



the observer to the source, and to work with Dh ^ 

dm) 

Dhy, instead. In most of the following, we convert 
from code units, which are in c = G = Mq = 1, to cgs 
units when stating and plotting numerical results. The 
conversion factors we use are 1 Mq = 1.4772 x 10^ cm for 
length, and 1 Mq = 0.004927 ms for time. For simplicity, 
we state the radii of GW extraction spheres and world- 
tube radii in code units. These and their corresponding 



cgs values are listed in Table III 



IV. RESULTS 

In this section, we compare the most reliable extrac- 
tion method that contains the least amount of systematic 
errors, CCE, with the various other curvature-based ex- 
traction methods, i.e., RWZM and NP extraction (both 
at finite radii). We also perform a comparison with the 
quadrupole approximation which has been employed in 
virtually all core collapse simulations to date. 

This section is structured as follows. First, we review 
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briefly the morphology of the gravitational waveforms ex- 
pected from rotating core collapse and bounce. 

Second, we elaborate in more detail on the method 
with which we obtain the gravitational strain h from the 
quantities measured during the simulation. This is im- 
portant, since the derived strain typically contains severe 
non-linear drifts making a proper analysis largely impos- 
sible without significant preprocessing. 

Third, we assess the accuracy of each individual 
method, i.e., we analyze the radial dependence of the NP 
and RWZM extraction methods, since they are strictly 
valid only in an asymptotic frame at an infinite distance 
from the source where any contributions from the stress- 
energy tensor vanish. Since the matter densities are non- 
zero at the CCE world-tube locations, we also analyze the 
radial dependence of the waveforms extracted via CCE. 

Fourth, we compare the results obtained via NP and 
RWZM extraction, and the approximate QF with results 
obtained via CCE. 

Finally, we perform a convergence check on the com- 
puted waveforms by using a set of three different resolu- 
tions. 



that decreases in amplitude as the core approaches its 
final equilibrium. 

While the overall morphology of the GWs emitted by 
the three models is the same, there are subtle differ- 
ences that are worth commenting on. Models A1B3G3 
and A3B3G3 produce so-called type-I signals [SU] with 
a single pronounced major peak at core bounce. Since 
model A3B3G3 is more rapidly spinning, its inner core is 
more deformed, and hence produces a stronger GW sig- 
nal at core bounce than model A1B3G3. Model A1B3G5 
has a very small inner core at bounce and produces a 
type-Ill signal that is characterized by a much less pro- 
nounced bounce peak and generally low-amplitude GW 
emission. Note that type-II signals, characterized by mul- 
tiple wide and pronounced bounce peaks seen in early 
work [501 ll03H106j have been demonstrated to disap- 
pear in simulations using general relativity and a proper 
electron-capture treatment [H [S] . 



B. Computing the Strain 



Morphology of Rotating Core Collapse 
Waveforms 



The core collapse models considered in this work re- 
main nearly axisymmetric during collapse and emit GWs 
predominantly via the even-parity (^, m) = (2,0) spher- 
ical harmonic mode. This mode has a maximum on the 
equator and, hence, we plot all waveforms as seen by 
an observer in the equatorial plane. We write /i+,e = 



,20 V 
1+ -2-'20 



: 7r/2,(/)) where _2^20 is the spin-weighted 
spherical harmonic with spin s = —2. Note that the 
{£, m) = (2, 0) mode is axisymmetric and thus, the equa- 
torial strain h^ ,, has no 0-dependence. 

We convert to cgs units by using the transformation 
as discussed in Sec. IIIE[ and we align the maxima of 
the waveforms such that they occur at i = 0.0 ms, cor- 
responding roughly to the time of core bounce in each 
model. 

The waveforms of the three models are shown in Figs.[l] 
(model A1B3G3), [2] (model A1B3G5), and |3] (model 
A3B3G3). All models exhibit a very similar behavior. 
Prior to core bounce {t < ms), the matter undergoes 
an aspherical accelerated collapse. Due to this aspher- 
ical acceleration, the GW signal is monotonically rising 
until it peaks when the contracting inner core is drasti- 
cally decelerated. This deceleration is caused by the sud- 
den stiffening of the EOS as a result of nuclear repulsive 
forces which emerge when nuclear densities are reached. 
During this deceleration, the GW signal becomes rapidly 
negative, reaching its second peak (the "bounce peak") 
roughly when the core rebounds. Subsequently, the inner 
core undergoes a relaxation phase (ring-down) in which 
it loses its remaining pulsation energy by launching sec- 
ondary shocks. This results in an oscillatory GW signal 



We first consider the computation of the strain /i+^e 
from the RWZM formalism. Since our models emit GWs 
predominantly in the {i,m) — (2,0) even mode , the 
computation of the strain from the even- and odd-parity 
RWZM master functions reduces to Eq. (50) so that no 
time integral is necessary to obtain h^ ,,. However, we 
still notice an unphysical drift in the waves. Since the 
RWZM master functions are computed at a finite dis- 
tance from the source, we have the following systematic 
errors (summarized as the "finite- radius error"): (i) a 
non- vanishing matter density at the extraction spheres, 
(ii) near-zone effects and (iii) gauge ambiguities. The lat- 
ter error arises as a result of deviations from the Bondi 
gauge (see [90| for an improvement). The artificial drift 
is part of the finite-radius error, since it is becoming less 
pronounced with increasing extraction radius. 

In order to reduce the contribution of these artifi- 
cial low-frequency components, we first transform to the 
Fourier domain, multiply by / in order to take the first 
time derivative, and then apply fixed FFI [102] to obtain 
/i+^e- This procedure effectively acts as a filter that sup- 
presses unwanted low-frequency components and at the 
same time minimizes spurious oscillations in the time do- 
main such as Gibbs ringing or additional non-linear low- 
frequency drifts. 



The QF (cf. Eq. 36) directly computes the strain and 



does not suffer from low-frequency drifts. However, the 
NP and CCE methods compute the second time deriva- 
tive of the strain and, hence, must be integrated twice 
in time (Eq. 66). For this, we employ FFI to minimize 
the influence of artificial low-frequency components. Un- 
fortunately, the time integration is still subject to some 
amount of low-frequency error as we shall discuss in the 
following. 
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1. Error in the Time Integration 

FFI introduces a free parameter /o that must be chosen 
based on the expected lowest physical frequency compo- 
nent of the signal. It must be larger than the spurious 
artificial low-frequency contributions that are introduced 
by aliased numerical noise and spectral leakage |102) . 

Unfortunately, in the case considered here, the arti- 
ficial low-frequency contributions overlap with the low- 
frequency part of the physical signal. Since it is not pos- 
sible to disentangle physical from artificial contributions 
at a given frequency, we have to choose an /o that is 
larger than the highest unphysical frequency contained 
in the signal. Thus, part of the overlapping physical low- 
frequency spectrum is lost when constructing the strain 
h. 

We identify the highest unphysical frequency by choos- 
ing a set of different /o for a given waveform, i.e., we 
introduce a family of strains h{t; /o), and by imposing a 
relative maximum deviation maxj (S/i(i; fo)/5fo between 
the resulting strains h(t; /q) during ring-down that is not 
larger than some small e. Since the ring-down phase is at 
the end of the wave train, the impact of the accumulated 
drift is largest here and can be clearly identified. In- 
creasing /o reduces the drift in the ring-down phase, but 
also removes physical content, i.e., the monotonic rise of 
the signal in the prebounce phase. In order to gauge 
how much information is lost prior to core bounce, we 
compute the differences of h{t; Jq) from the quadrupole 
waveform in an interval t G [—10 ms, ms], since the 
quadrupole waveforms do not suffer from time integra- 
tion issues and are presumably accurate up to at least 
the late prebounce phase. If we stop at some level of tol- 
erance for any deviations maxt Sh(t; /o)/^/o during ring- 
down and deviations from the quadrupole waveform prior 
to core bounce, we have identified an appropriate /q. In 
practice, we choose a threshold maxt6Dh{t; /o)/^/o ^ 
e ~ 0.1 cm/Hz during ring-down. 

Our numerical experiments show that the cut-off fre- 
quency /o is model and extraction-method dependent 
and must be determined individually for any new set of 
initial data. In Table [III li^t the frequencies /o for each 
of the considered models and extraction methods which 
yield the lowest deviations during ring-down and at the 
same time resemble as closely as possible the quadrupole 
waveform in the prebounce phase. 

We find that NP and RWZM extraction, which both 
operate at finite radii, are subject to stronger drifts than 
the CCE method which computes the waveforms gauge 
invariantly at future null infinity J^'^. This is not sur- 
prising, given that the two former methods both suffer 
from near-zone and gauge errors which typically lead to 
low-frequency drifts in the waveform. Hence, the strain 
h'^^ as computed by CCE retains most of the physi- 
cal information at the low-frequency end of the spectrum 
with a cut-off at /o = 100 Hz. Unfortunately, even this 
value may not be low enough, given that this frequency 
falls right into the band of highest sensitivity of km-scale 



TABLE II: Lowest possible physical frequencies that result in 
strain amplitudes with deviations maxt \SDh(t; fo)/Sfo\ of no 
more than e ~ 0.1 cm/Hz. In all cases, CCE yields the lowest 
possible /o and, hence, retains most physical information at 
the low-frequency end. 





A1B3G3 


A1B3G5 


A3B3G3 


Method 


/o [Hz] 


fo [Hz] 


fo [Hz] 


NP (i?cx = IOOOM0) 


300 


300 


250 


RWZM (iicx = IOOOMq) 


250 


400 


200 


CCE {Rr = IOOOM0) 


100 


100 


100 



ground-based detectors |1071 1108] . Not being able to re- 
solve the low-frequency components is clearly a drawback 
of the curvature-based extraction methods. 

Fortunately, as we will discuss in more detail in 
Sec. |IVD H frequencies below 100 Hz do not contribute 
significantly to the inferred theoretical signal-to-noise ra- 
tio (SNR) for the models considered in this study or for 
the GW signal associated with core bounce in slowly to 
moderately rapidly rotating core collapse in general [S]. 
Hence, at least the CCE method yields robust predictions 
for detection. A closer and more detailed comparison be- 
tween the waveforms computed with the various methods 



will be performed in Sec. IV D 1 



In the following, we use the cut-offs fo as given in Ta- 
ble |TT] for the various models and extraction methods. 



C. Radial Dependence 

The physical gravitational strain h scales oc 1/D with 
distance D, provided it is computed in an asymptotic 
frame at large distances from the source D 00 (i.e., at 
astrophysical distances). At large asymptotic distances, 
wc should therefore observe Dh — const. Since NP ex- 
traction and the RWZM formalism are evaluated at a 
finite distance from the source, they are both subject 
to finite-radius errors and we will generally not exactly 
observe Dh — const. Rather, we expect the signal to 
converge with increasing extraction radius towards its 
asymptotic shape and magnitude. In the context of vac- 
uum binary black hole mergers, the asymptotic behavior 
and finite-radius error of NP extraction has been ana- 
lyzed in Ref. [109], where it is found that extrapolations 
based on extractions at radii R > 300 M, where M is 
the total mass of the system, yield acceptable results. 
In our case, however, the finite-radius error contains the 
additional error arising from non-zero matter content at 
the CCE world-tube locations and NP/RWZM extrac- 
tion spheres. 

In Table |III[ we summarize the time-averaged matter 
densities (p) at the various extraction spheres in our mod- 
els. For simplicity, we do not compute the extraction- 
surface-averaged matter density but simply report the 
equatorial density at the various extraction spheres. The 
most compact model A3B3G3 has (p) that are a factor 
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of a few smaller at any given radius than in the other 
models. Therefore, we expect the systematic finite-radius 
error to be smallest in model A3B3G3. 

In order to quantify the finite-radius error, we com- 
pute Dh+^e at a succession of extraction spheres with in- 
creasing radii i?ex = {500 Mq , 700 Mq , 800 Mq , 900 M© , 
1000 Mq} and evaluate the differences. For R^x > 
IOOOM0, the spatial resolution of our computational grid 
becomes too coarse for accurate wave extraction and we 
limit our analysis to i?cx < IOOOMq (Table |III|. Note 
that for a given model and extraction method, we use 
the same cut-off frequency /o for all extraction radii and 
world-tube locations. 

In principle, we should extrapolate the waveforms as 
obtained at the different extraction spheres to infin- 
ity. We observe, however, that the differences at large 
radii are within our numerical errors (see Sec. IV E) and 
asymptote rapidly. Therefore, we simplify the analysis 
by inspecting the behavior at successive radii without 
extrapolating. We will see in Sec. IV Dl that this ap- 
proach is justified. The CCE method, which does not 
require any extrapolation, shows good agreement with 
results obtained at finite radius within our numerical er- 
rors. 



1. The NP Scalar ^4, 

In the upper left panels of Figs, [l] [2j and [3) we 

show Dh^^^ as computed from the NP scalar \I'4 for 
model A1B3G3, A1B3G5 and A3B3G3, respectively. In 
the bottom panel, we show the absolute differences a 
between -D/i+^ at Rex from the reference distance at 
Rex — 1000 Mq. In an ideal asymptotic frame, all curves 
would line up exactly. This is not the case in practice. 
We notice that the curves asymptote with increasing ex- 
traction radius, i.e., the differences a between two suc- 
cessive extraction spheres converge to zero. This behav- 
ior shows that our extraction radii, albeit rather close 
to (and even inside) the star, lead to finite-radius errors 
for the waveforms computed from the NP scalar ^'4 that 
are below the discretization errors (cf. Sec. IV E). Wc 
measure a maximum absolute difference in amplitude of 
(T < 2 cm between the two outermost extraction spheres 
of any model (see lower panels of Figs.[Tj[2j and[3|. This 
corresponds to a relative error of no more than ~ 2% 
when compared to the peak amplitudes. 

Note that we had to cut-off low frequencies (Table [lT|) 
in order to remove spurious non-linear drifts. This re- 
stricts our analysis to frequency components above the 
cut-off frequency /q. Experiments show that for larger 
extraction radii, artificial low-frequency components be- 
come less prominent. Hence, extracting at greater ra- 
dius would allow us to decrease /o, but is presently too 
computationally demanding to be possible for production 
simulations. 



2. The Regge- Wheeler- Zerilli-Monerief Formalism 



The RWZM variables are computed as perturbations 
on an assumed fixed background spacetime. In an asymp- 
totic frame at large distances from the source, they are 
independent of radius. As in the case for NP extrac- 
tion, this is not achieved in practice, but the residual 
errors should converge with increasing extraction radius. 
We measure a relative difference in amplitude between 
the two outermost detector spheres of ct < 17 cm for all 
models (see the right panels of Figs. [11 [2] and |3| . This 
corresponds to a relative error of < 8% when compared 
to the maximum amplitudes and is significantly larger 
than what we find for NP extraction. 

In addition, the RWZM method produces high- 
frequency variations in the waveform at core bounce and 
similar high-frequency features in the ringdown phase 
that are not seen in GWs extracted with the other meth- 
ods. These features do not appear to converge with in- 
creasing radius; at least not at radii accessible to our sim- 
ulations. They are particularly manifest in GWs of mod- 
els producing weak signals, e.g., in the signal emitted by 
model A1B3G5 of our model set (Fig. [2]). Furthermore, 
and most pronounced in model AlB3G3's waveform, a 
large spike during core bounce is visible in the RWZM 
result, but is not produced by any of the other meth- 
ods (see Fig. [l] and the comparison in the upper panel of 
Fig.|4|. 

In order to investigate the cause of the differences seen 
with RWZM extraction, we perform a range of test cal- 
culations. These include (j) using two additional in- 
dependent implementations of the RWZM method, one 
assuming a Minkowski background, the other using a 
generalization of the RWZM approach |110j . (iz) per- 
forming a computationally very expensive simulation 
with extended grids, allowing for RWZM extraction at 
Rex — 3000 Mq, (m) performing simulations with up to 
a factor of 2 higher resolution and modified mesh refine- 
ment boundary locations, and {iv) changing the space- 
time gauge conditions, including exponential damping of 
the evolution of the coordinate shift at large radii near 
the extraction spheres. 

None of the above tests leads to any significant change 
of the RWZM result. This brings us to the conclu- 
sion that the high-frequency features observed in RWZM 
waveforms are systematic problems tied, most likely, to 
the particular perturbative nature of the RWZM scheme. 
One notable difference of the RWZM formalism from the 
other methods is the procedure of projecting out the 
spherical background geometry (e.g. Eq. 48 1. This can 



result in very small values for the aspherical perturba- 
tion coefficients that are prone to numerical noise and 
cancellation effects. The RWZM approach may therefore 
be less suitable for the extraction of the generally weak 
GW signals emitted in core collapse. 
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TABLE III: Time-averaged equatorial matter densities and their variations at extraction radii for the computed models. The 
first two columns report the extraction radii in code units A/q and in cgs units, respectively. All extraction surfaces are 
located on the fourth refinement level with a spatial resolution of Ax = 8M0 (~ 11.82 km) and time resolution of Ai = I.6M0 
(~ 7.9 X 10-*^ s). 



Rex 




(p) A1B3G3 


(p) A1B3G5 


{p) A3B3G3 




[km] 


[g cm-3] 


[g cm-3] 


[g cm-3] 


500 


739 


(1.2 ±0.3) X 10** 


(1.6 ±0.4) X 10** 


(8.6 ±0.2) X 10** 


700 


1034 


(2.3 ±0.4) X lO'^ 


(2.6 ±0.5) X 10'' 


(1.3 ±0.3) X 10'' 


800 


1182 


(8.6 ± 1.2) X 10** 


(9.3 ± 1.2) X 10'^ 


(3.5 ±0.9) X 10*^ 


900 


1329 


(3.0 ±0.4) X 10** 


(3.3 ± 0.4) X 10*^ 


(5.5 ± 1.6) X 10^ 


1000 


1477 


(9.3 ± 0.6) X lO'"^ 


(1.2 ±0.6) X lO'^ 


(2.5 ±2.2) X 10^ 




t [ms] t [ms] t [ms] 

FIG. 1: Top panel (from left to right): D/i+^e computed using NP extraction at radh _Rcx = (500, 700, 800, 900, 1000) M©, 
Dhl^^ computed with CCE at J+ using world-tube data at Rr = (500, 700, 900, 1000) Mq, and Dh^J'^ computed using 
the RWZM formahsm at radii Tiex = (500, 700, 800, 900, 1000) Mq, all for model A1B3G3. Bottom panels: Absolute difference 
a of the waves extracted at the various extraction radii and world-tube locations from those extracted at the outermost 
radius/location. The waveforms converge with increasing extraction radius and world-tube location. For NP extraction, we 
measure at the outermost detector sphere 7?ex = 1000 Mq a maximum difference to the next closest detector 7?ex = 900 Mq 
of (T4 = 2 cm, which corresponds to a percentage error of ~ 1% relative to the maximum. For CCE, we measure a maximum 
difference of 0-4 = 1.5 cm, corresponding to a percentage error of ~ 0.7% relative to the maximum. For RWZM extraction, we 
have (T4 = 14 cm, corresponding to ~ 4% relative to the maximum. We note that RWZM is subject to additional high-frequency 
features and also contains a spurious spike in the bounce peak. 



3. Radial Dependence on World- Tube Location for CCE 

The CCE method uses metric data from a time succes- 
sion of finite-radius coordinate spheres, the world-tube, 
as inner boundary data for the evolution of the gravita- 
tional field out to . In vacuum, the method does not 
depend on the particular choice of any given world-tube 
radius [371 HH]- However, the presence of matter at the 
world-tube location leads to a systematic error and im- 
poses an artificial dependence of the waveforms computed 
at on the world-tube location. We plot in the center 
two panels of Figs, [l] [2j and [3] the waveforms obtained 



at J'^ using different world-tube radii as inner bound- 
aries for the characteristic evolution. The center bottom 
panel of these figures depicts the absolute difference a 
between the waveforms from the outermost world-tube 
radius at i?r — 1000 Mq and the waves from each of 
the smaller world-tube radii. It is apparent that the dif- 
ferences between the outermost two world-tube radii is 
always smallest, with absolute differences a < 1.6 cm for 
all models and all times, and with an error relative to 
the maximum amplitude of < 1% for models A1B3G3 
and A3B3G3, and - 6% for model A1B3G5. For the lat- 
ter model, we also notice strong drifts at the innermost 
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FIG. 2: Radial dependence of waveforms computed for model A1B3G5. See caption of Fig. [T] for details. For NP extraction, 
we measure at the outermost detector sphere _Rcx ~ 1000 M© a maximum difference of a4 = 0.4 cm to the next closest detector 
Jlex = 900 M0. This corresponds to a percentage error of ~ 2% relative to the maximum. For CCE, we measure a maximum 
difference of 0-4 = 1.6 cm, corresponding to a percentage error of ~ 6% relative to the maximum. For RWZM extraction, we 
have 0-4 = 2.4 cm, corresponding to ~ 8% relative to the maximum. We note that RWZM is subject to strong additional high- 
frequency features. We also note that, while CCE shows for this model a larger error due to matter effects than finite-radius 
'I'4 extraction, it is more accurate at low frequencies (cf. Table [ll| . 




-2 2 4 
t [ms] 



-2 2 4 
t [ms] 



10-4 
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FIG. 3: Radial dependence of waveforms computed for model A3B3G3. See caption of Fig. [T] for details. For NP extraction, 
we measure at the outermost detector sphere iicx = 1000 Mq a maximum difference of 0-4 = 2 cm to the next closest detector 
7?ox = 900 Mq. This corresponds to a percentage error of ~ 1% relative to the maximum. For CCE, we measure a maximum 
difference of a4 = 1 cm, corresponding to a percentage error of ~ 0.2% relative to the maximum. For RWZM extraction, we 
have (74 = 17 cm, corresponding to ~ 4% relative to the maximum. 
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world-tube radii. Systematic errors due to the presence 
of matter can therefore become significant when the sig- 
nal is weak and the density large. Note that the strong 
drift at the innermost world-tubes may be removed by 
an increased FFI cut-off frequency of /o = 150 Hz in this 
case. Since we would like to retain as much physical 
information as possible, and since the outermost world- 
tube location permits /o = 100 Hz, we have chosen this 
value for all world-tubes. Generally, the lower FFI cut- 
off frequency /o induces a more sensitive radial behavior 
with respect to low-frequency drifts. Since CCE permits 
a lower /o than NP extraction, the radial variations are 
slightly larger for CCE in model A1B3G5 when compared 
to the radial variations of NP extraction (cf. bottom pan- 
els of Fig. [2|. In the other models, we find smaller radial 
variations between the two outermost CCE world-tubes 
than between the two outermost NP extraction spheres 
(cf. bottom panels of Figs, [l] andjs]), even though fo is 
smaller for CCE than for NP extraction in these mod- 
els as well. This indicates that the remaining systematic 
non-zero matter error in CCE is not as important as the 
additional near-zone errors and gauge ambiguities inher- 
ent to NP extraction. 



D. Comparison 

1. Comparison with Cauchy- Characteristic Extraction 

CCE yields waveforms that contain the least amount 
of systematic errors compared to the other extraction 
methods considered in this work. It is completely gauge 
invariant and is free of near-zone errors. As found in 
Sec. |IVB 1[ it is the only curvature-based method that 
captures most of the low-frequency band. Furthermore, 
as found in Sec. |IV C[ the remaining error due to the 
non-zero stress-energy tensor is small. We therefore use 
the waveforms obtained with CCE as a benchmark. In 
Figs. |4j [5]and[6j we examine the waveforms as obtained 
by the various extraction methods, i.e., NP, RWZM, QF, 
and CCE for each model. In each figure, the top panel 
displays the amplitudes of the waves, D/i+^e, with the 
panel right below showing the absolute differences a of 
each extraction method from CCE. It is apparent that in 
all cases the RWZM formalism yields the largest differ- 
ences from CCE. As discussed in Sec. |rVC| the RWZM 
formalism is subject to high-frequency noise and yields a 
large spurious spike at core bounce, most pronounced in 
model A1B3G3. 

The quadrupole approximation and NP extraction 
only lead to small differences to CCE, in particular at 
core bounce, but also in the ring-down part of the wave- 
form. Moreover, it is surprising that the quadrupole for- 
malism yields agreement so remarkably close to the re- 
sults obtained via CCE, given the rather simplistic na- 
ture of the QF. Quantitatively, when compared to the 
waveforms obtained via CCE, we find for model A1B3G3 
that the waves obtained via the quadrupole approxima- 
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FIG. 4: Comparison of waveform amplitudes Dh+^e, their ab- 
solute differences a from CCE waveforms, and spectrograms 
of the power ratio Lds between waveforms obtained from 
each extraction method and waveforms obtained via CCE 
for model A1B3G3. If LdB = 0, the waveform of the par- 
ticular extraction method yields equal power per time and 
frequency bin with respect to that obtained with CCE. This 
is indicated by green colors. Blue colors indicate less power, 
and red colors indicate more power. See text for details. NP 
extraction and the quadrupole approximation yield remark- 
able agreement with CCE at frequencies below 2 kHz at and 
after core bounce, while the RWZM formalism yields a spuri- 
ous spike during core bounce and generally contains artificial 
high-frequency oscillations. This also becomes clear in the 
spectrograms of the power ratio between RWZM formalism 
and CCE since LdB > over a wide range of time and frequen- 
cies {bottom panel). Prior to core bounce —8ms < t < —1ms, 
NP extraction results in less power compared to CCE, while 
the QF yields more power. 



tion lead to smaller core-bounce peaks with differences 
up to ^ 10cm 5%). For the same model, NP extrac- 
tion results in core-bounce peaks that are smaller as well, 
with differences of ~ 5%. Note, however, that in the NP 
waveform, the first positive peak prior to bounce is much 
larger ~ 31% than what is predicted by CCE. The QF 
result, on the other hand, agrees much better with CCE 
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FIG. 5: Comparison of waveform amplitudes Dh+,e, their ab- 
solute differences a from CCE waveforms, and spectrograms 
of the power ratio Lds between waveforms obtained from 
each extraction method and waveforms obtained via CCE for 
model A1B3G5. If L^b = 0, the waveform of the particular 
extraction method yields equal power per time and frequency 
bin with respect to that obtained with CCE. This is indicated 
by green colors. Blue colors indicate less power, and red col- 
ors indicate more power. See text for details. The waveforms 
from NP extraction and quadrupole approximation agree well 
at frequencies below 2 kHz at and after core bounce, while the 
RWZM formalism is subject to artificial high-frequency oscil- 
lations. The spectrograms of the power ratio between wave- 
forms from the RWZM formalism and CCE in the bottom 
panel further support this, since Lds > in this case. Prior 
to core bounce —8 ms < t < —1ms, NP extraction results 
in less power compared to CCE, while the QF yields more 
power. 



at this peak (~ 5% overprediction) . 

For model A3B3G3, we find that at the bounce peak 
the QF (NP) amphtudes are 11% 7%) smaller than 
the CCE prediction. For this model, the first positive 
peak before bounce is overpredicted by NP by ^ 69%, 
while the QF yields an overprediction of only ^ 11% 
compared to CCE. 

A separate treatment is necessary for model A1B3G5. 
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FIG. 6: Comparison of waveform amplitudes Dh+^e, their ab- 
solute differences a from CCE waveforms, and spectrograms 
of the power ratio Lds between waveforms obtained from 
each extraction method and waveforms obtained via CCE 
for model A3B3G3. If LdB = 0, the waveform of the par- 
ticular extraction method yields equal power per time and 
frequency bin with respect to that obtained with CCE. This 
is indicated by green colors. Blue colors indicate less power, 
and red colors indicate more power. See text for further ex- 
planation. The waveforms from NP extraction and from the 
quadrupole approximation agree well below 2 kHz, while the 
waveform from the RWZM formalism contains artificial high- 
frequency oscillations thus leading to a higher power emission 
than waveforms obtained via CCE {bottom panel). Prior to 
core bounce —10 ms < t < —2ms, NP extraction results in 
less power compared to CCE, while the QF yields more power. 



As briefly discussed in Sec. |IVB 1[ the physical low- 
frequency components are filtered out by FFI in 
curvature-based extraction methods. The waveform of 
model A1B3G5 is most affected by this: Any physi- 
cal low-frequency modulations or offsets are removed. 
Hence, the waveforms from curvature-based extraction 
are shifted downwards with respect to the QF waveform 
that does not require filtering. This shift leads to large 
absolute differences in the peak amplitudes by as much 
as ~ —10 cm at the first and second peaks when com- 
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pared to the QF waveform, yielding relative differences 
by as much as ~ 90% at core bounce. Constant (or 
nearly constant) offsets in the waveforms, however, are 
not visible to GW detectors. In order to get a better 
measure of the differences between the various extraction 
methods, we therefore compare the change in amplitude 
SDh+^f. = \Dh*^ g — -D/i^ g| between first and second peak 
occurring at time ti and t2, respectively. Compared to 
CCE, the change 5Dh!^^ measured in the quadrupole 
waves is ^ 7% smaller and 5DhJ^^ of the NP waveform 
is ~ 1% larger. 

Since a hypothetical matched-filtering GW search 
would be sensitive to differences in phase, an important 
aspect is the phase relation of the waveforms. A measure 
for the phase is given by the time lag between succes- 
sive wave peaks, and we label the time of the peaks by 

= ii, ^2, ^3, • ■ • in their temporal order of occurrence. 
Note that in all cases, we have aligned the waveforms at 
the second peak occurring at t2 = ms so that the time 
lag at this peak is zero for all methods and all models. We 
therefore measure the time lag 5tn^2 for a peak i„^2 rela- 
tive to the second peak, e.g., we measure 5ti^2 = 1^2 —ti\- 
By comparing a number of time lags i5fn,2j we gener- 
ally find that NP extraction, CCE, and QF produce the 
same phasing with an error of less than ^ ±0.05 ms for 
all peaks and all models. The differences in the time lags 
between successive wave peaks computed from the var- 
ious extraction methods are therefore close to our time 
resolution of Ai « 8 x 10"^ ms. Note that the RWZM 
formalism is excluded from this analysis, since the addi- 
tional high-frequency components make it hard to clearly 
identify the times of the maxima. Visual inspection sug- 
gests a phase error for RWZM comparable to the other 
methods, provided high-frequency contributions are ig- 
nored. 

In the bottom three panels of Figs. |4j [5] and |6j 
we plot spectrograms'^ of the power ratios LdB — 
101ogiQ(Fi/Fo)- Here, the power spectrum Pi = \h{f)\^ 
is computed from the quadrupole approximation, NP ex- 
traction, and the RWZM formalism, respectively. The 
power spectrum Pq — |^cce(/)P is computed from the 
waves obtained via CCE. At a given time i, the power 
spectra Po,i a-re obtained from the short-time Fourier 
transform of the strain over a time window of 2 ms cen- 
tered at t. Thus, LdB measures the power ratio per time 
and frequency bin of all extraction methods relative to 
the CCE method. If LdB = 0, the extraction method 
emits equal power per time and frequency bin and, hence, 
is equivalent to the waves from CCE. This is indicated 
by green colors. Red regions indicate that waves ob- 
tained with the corresponding extraction method emit 
more power than the CCE waves during that time and 
frequency bin; blue indicates less power. 



By inspection of the spectrograms, it is apparent that 
at core bounce, the NP method, and also the quadrupole 
approximation predict waves carrying roughly equal pow- 
ers with respect to the waveforms obtained via CCE in 
a time interval of [—4 ms, 4 ms] and over the entire fre- 
quency range. Furthermore, we observe that in the pre- 
bounce phase [—10 ms, 4 ms] , the quadrupole waves emits 
more power in low frequencies than the CCE waves, and 
the NP extracted waveforms emit less power. This is 
largely an effect of the different cut-off frequencies intro- 
duced in the integration of the waves obtained from the 
curvature-based methods (Table |ll|. 

We also find from the spectrograms that the RWZM 
extraction always deviates strongest from all other meth- 
ods, primarily because its GWs contain spurious addi- 
tional high-frequency components (cf. bottom panels of 
Figs. [4] [5 and [6]). This is clearly visible in the spectro- 
grams, which show additional red "speckles" of higher 
emitted power throughout the wavetrain and frequency 
band. 

As opposed to the quadrupole approximation, all 
curvature-based extraction methods are filtered below a 



frequency /o (Sec. IVB 1 1. It is crucial to know whether 
the missing low frequencies can spoil the detectability 
of the GWs. To gauge the influence of low-frequencies 
on the theoretical signal strength in GW detectors, we 
compute the theoretical optimal SNR for the quadrupole 
waveform once including low-frequency contributions and 
once artificially cutting them off. The SNR is given by 

m 



(69) 



^ The spectrograms are made up of 100 time bins of 0.2 ms each 
and use a Hann window with a width of 2 ms. 



where h is the Fourier transform of the strain h as mea- 
sured at a distance D = lOkpc to the source, and Sh 
is the one-sided noise power spectral density for a given 
detector, i.e., the detector sensitivity function. 

Table |IV| lists for all models the theoretical optimal 
SNRs of waves extracted with the quadrupole formal- 
ism for the LIGO |107l 1112] detector and for the zero- 
detuning high-power configuration of advanced LIGO 
[TU5] . By cutting off at /o = 100 Hz ("QF, /o = 100 
Hz"), we observe no significant reduction in SNR for 
models A1B3G3 and A3B3G3 and, hence, no loss of de- 
tectable GW information in rotating core collapse for 
these models. Model AlB3G5's type-Ill waveform, on 
the other hand, has low-frequency components of high 
relative strength. In this model, a cut-off at /o = 100 Hz 
already leads to a SNR reduction by ~ 17% (~ 23%) in 
LIGO (advanced LIGO). 

When using 100 Hz < fp < 300 Hz as a cut-off, we find 
a reduction in SNR as large as ^ 10 — 30% for all models. 
Comparing to the cut-off frequencies required for FFI 
( Table [ll|, this indicates that the finite-radius curvature- 
based NP and RWZM methods suffer from their inability 
to properly resolve frequencies in the range 100 Hz < f < 
300 Hz. 
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As shown in Table [TV] the theoretical optimal SNRs as 
computed from the waves obtained for models A1B3G3 
and A3B3G3 via CCE are essentially unaffected by low- 
frequency removal. Instead, they yield a slightly in- 
creased SNR compared to the QF result. This is due 
mostly to the higher amplitudes of the CCE waves at core 
bounce. Since low-frequency components are missing, 
model A1B3G5 is subject to a loss in SNR. CCE yields 
SNRs for LIGO and advanced LIGO that are identical 
to those obtained with QF with a cut-off at /o — 100 Hz 
("QF, /o^lGOHz"). 

The RWZM method always overpredicts the SNR by 
~ 40-100%. The reason is apparent from the spectro- 
grams in the bottom panels of Figs. |4j [5] and [6) in which 
spurious additional high-frequency components lead to 
artificially high GW power and a corresponding overes- 
timate of the SNR. 

The waves computed with NP extraction show a high 
degree of agreement in SNR with the waves obtained with 
the QF. This, however, is misleading, since there are two 
balancing effects: (i) NP extraction predicts a higher am- 
plitude in the first peak prior to core bounce, which would 
yield a larger SNR, but (ii), in NP extraction we must 
cut off frequencies / < 300 Hz, which artificially reduces 
the SNR. 

Finally, we address the question of whether the waves 
obtained with the various extraction methods are within 
the tolerance for detection in a (hypothetical) matched- 
filtering GW data analysis of LIGO and advanced LIGO. 
Ideally, waveforms for the same model lead to a detection 
of the same model parameters and hence should not be 
distinguishable within a given threshold. 

The (dis) agreement of waveforms obtained from dif- 
ferent methods can be quantified by the mismatch (see, 

e.g. [niiiii]) 

M,„is = 1 - M , (70) 
where the hesi match M. is given by 

= max max max ©[/ii, ^12] , (71) 

*0 01 02 

which involves a maximization over time of arrival and 
the two phases 0i and (^2 of the two wave signals h\ and 
/i2j respectively. The overlap O between two waveforms 
is given by 



{hi\h2 



^{h,\h,){h2\h2) ' 

with the detector-noise weighted scalar product 

Mmif) 



(/ii|/i2) =4Re / df- 
Jo 



Shif) 



(72) 



(73) 



hi 



A mismatch jMmis of zero indicates that waveforms 
and /i2 are identical. Conversely, a mismatch of Mmis — 
1 indicates that the waveforms are completely different. 
In Table [V] we list the mismatches between each of the 



TABLE IV: SNR ps according to Eq. (|69| for all models and 
all extraction methods at a distance D = lOkpc for LIGO 
(top) and advanced LIGO (bottom). Note that we take into 
account only those frequencies above a given /o that have been 
determined to be reliable (see Table [ll|. Since the quadrupole 
formalism is robust at low frequencies, we can gauge the influ- 
ence of the neglected frequencies by additionally computing 
the SNR for the quadrupole waveform with the same low- 
frequency cut-off. For LIGO, we find no loss in SNR. [] For 
advanced LIGO the loss is at ~ 1%. An exception is model 
A1B3G5 where the low-frequency components contribute sig- 
nificantly to the total emission. 



Method 


A1B3G3 


A1B3G5 


A3B3G3 




ps 


PS 


ps 


LIGO 








QF 


3.5 


0.6 


7.8 


QF (/o = 100 Hz) 


3.5 


0.5 


7.8 


NP (i?cx = lOOOA'/o) 


3.1 


0.4 


7.4 


RWZM (iiex = IOOOMq) 


6.9 


0.7 


17.0 


CCE {Rr = IOOOMq) 


3.8 


0.5 


8.5 


advanced LIGO 








QF 


49 


9 


95 


QF (/o = 100 Hz) 


49 


7 


94 


NP (i?cx = IOOOM0) 


49 


6 


96 


RWZM (7?ex = IOOOMq) 


105 


12 


209 


CCE {Rr = IOOOMq 


52 


7 


103 



TABLE V: Mismatch Mmia according to Eq. ( 70 1 for all mod- 
els between CCE and the other extraction methods for the 
LIGO (top) and advanced LIGO detector (bottom). Note 
that we take into account only those frequencies above a given 
/o that have been determined to be reliable (see Table |lT]| . 
The waveforms from the quadrupole approximation yield the 
smallest mismatch to the waves from CCE, except in model 
A1B3G3. But note that the quadrupole approximation yields 
waveforms that allow the inclusion of lower frequencies than 
all other methods and hence allow the computation of the 
mismatch over a greater frequency range (only limited by the 
cut-off frequency of CCE) which introduces a small bias. 



Method 


A1B3G3 


A1B3G5 


A3B3G3 










LIGO 








NP (i?cx = IOOOMq) 


5 X 10"^ 


19 X 10"^ 


8 X 10"^ 


RWZM (iJex = IOOOMq) 


12 X 10"^ 


38 X 10"^ 


5 X 10"^ 


QF 


6 X 10"^ 


13 X 10"^ 


2 X 10'^ 


advanced LIGO 








NP (i?ox = IOOOMq) 


3 X 10"^ 


12 X 10"^ 


7 X 10"^ 


RWZM (iiex = IOOOMq) 


14 X 10"^ 


48 X 10"^ 


7 X 10"^ 


QF 


4 X 10"^ 


7 X 10"^ 


2 X 10"^ 



extraction methods and the CCE method for all mod- 
els. Note that we compute the mismatch starting from 

fo = max{/Q"^\ /q^^}, where /q^^ and /g^^ are the lower 
cut-off frequencies as listed in Table [H] for the waveforms 
hi and /12, respectively, since we do not trust waveforms 
below their value of fo- We find that in all cases, the 
quadrupole approximation agrees best with waveforms 
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obtained via CCE with mismatches to within 1% or bet- 
ter: We find a mismatch in the quadrupole waveforms 
to CCE of 0.6% (0.4%) for LIGO (advanced LIGO) for 
model A1B3G3, 1.3% (0.7%) for LIGO (advanced LIGO) 
for model A1B3G5, and 0.2% (0.2%) for LIGO (advanced 
LIGO) for model A3B3G3. As reported in Table |v) NP 
extraction leads to slightly larger mismatches, due mainly 
(i) to less emitted power in the low-frequency band at and 
above this method's cut-oflF frequency, and (ii) to higher 
emitted power in the first wave peak (cf. spectrograms in 
second-lower panels of Figs. [Ij [sj and |6f . The RWZM 
formalism performs worst for models producing weak 
signals, e.g., model A1B3G3 and, in particular, model 
A1B3G5. This is due primarily to the artificial high- 
frequency components produced by this method. Note 
that the mismatch discussed above depends on the cut- 
off frequency /q. As a result, the range of frequencies 
contributing to the mismatch calculation is greatest for 
CCE-QF and smallest for CCE-RWZM. Hence, a full 
unbiased one-to-one comparison of the computed mis- 
matches is not possible. 

We next investigate the implications of the mismatch 
on detecting a particular model in matched-filtering anal- 
ysis. A reduction in the match M is equivalent to a re- 
duction in strain amplitude of the exact signal h (which 
here we assume to be given by the waveform computed 
via CCE) by A4h, hence effectively reducing the range 
of a GW detector by a factor of A^. To zeroth-order 
approximation, the number of detected events is propor- 
tional to the range cubed. A reduction in range by A4 
means a reduction of the number of detectable events by 
Ai^. If we require a loss of no more than 10% of all de- 
tectable events, the match (mismatch) between template 
waveform and exact signal must therefore never go be- 
low (above) M = 0.965 (TW^is = 3.5 x 10"^) pl3UTT5] . 
This indicates that when used as hypothetical templates 
in matched-filtering analysis of the LIGO and advanced 
LIGO data stream, NP extraction and the quadrupole 
approximation yield waveforms that are within the error 
tolerance, but RWZM is generally not. 

Overall, we conclude that NP extraction performs 
slightly worse than the quadrupole approximation when 
compared to CCE. The main reasons for this are (i) that 
NP requires a higher low-frequency cut-off and is there- 
fore missing important low-frequency components, (ii) 
that NP yields larger values for the first wave peak com- 
pared to what is obtained with CCE or the quadrupole 
approximation, and (iii) that the mismatches between 
NP and CCE waveforms are larger than those between 
the QF and CCE waveforms. The RWZM formalism is 
generally performing the worst since it produces artificial 
high-frequency contributions. 



2. Variations of the Quadrupole Formula 

In Figs. [7j |8] and [9] we plot waveforms for models 
A1B3G3, A1B3G5 and A3B3G3, respectively, computed 
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FIG. 7: Comparison of quadrupole to CCE waveforms for 
model A1B3G3. At and immediately after core bounce, the 
PV variant leads to a marginally smaller absolute difference 
to CCE (~ 1%) than any of the other QF variants. 



via t he QF given in Eq. (36 1, its PV and VS variants 
(Sec. Ill A I, and the waves as predicted by the CCE 



method. In the lower panel of these figures, we plot the 
absolute differences a from the CCE method for each 
variant of the QF. At core bounce, the smallest difl^er- 
ence from CCE for model A1B3G3 is predicted by the 
PV variant (^ 1% overprediction) , followed by the stan- 
dard QF (^ 5% underprediction) , and finally the VS 
variant (~ 8% underprediction). For model A3B3G3, 
we measure the smallest difference from CCE in the PV 
variant (~ 5% underprediction), followed by the stan- 
dard QF (~ 11% underprediction), and the VS variant 
(~ 13% underprediction). For model A1B3G5, the cut- 
off of low-frequency components leads to an offset of the 
CCE waveform compared to the waves obtained with the 
QF variants. We therefore compare the change in ampli- 
tude SDhj^ f, between the first and second peaks. We find 
that when compared to CCE, the change SDh^^ is ^ 2% 
smaller in the PV variant, ~ 7% smaller in the standard 
QF, and also ~ 7% smaller in the VS variant. Overall, 
the waves computed with the PV variant are closest to 
the results obtained with the CCE method. Since the 
definition of the QF is ambiguous, this finding may de- 
pend on the particular system studied and we cannot 
make strong general statements in support of one or the 
other variant. 
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FIG. 8: Comparison of quadrupole to CCE waveforms for 
model A1B3G5. The PV variant results in the smallest dif- 
ferences of the change SDh^^^ between the first and second 
peaks (by ~ 2%) compared to CCE. 
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FIG. 9: Comparison of quadrupole to CCE waveforms for 
model A3B3G3. As for model A1B3G3 (Fig.[7|, the PV vari- 
ant leads to the smallest difference from CCE (~ 5%). 



E. Convergence 

We check for convergence of our results via a resolu- 
tion study of model A3B3G3 using three different res- 
olutions, with finest resolutions Ax — 0.3 Mq (low), 
Aa; = 0.25 Mq (medium; our baseline resolution), and 
Aa; = 0.20 M© (high). 

In (relativistic) hydrodynamics simulations, conver- 
gence is notoriously difficult to analyze. The reasons are 
two-fold: First, the occurrence of hydrodynamical shocks 
reduces the accuracy locally to first order. In our models, 
a shock appears right after core bounce and has signifi- 
cant impact on the order of convergence of our scheme. 
Second, our simulations are subject to some turbulence 
of the fluid motion appearing soon after bounce. Turbu- 
lence is stochastic in nature. Even a slight change of the 
resolution can result in completely different local behav- 
ior of the fluid. For this reason, it is impossible to check 
convergence locally at each grid point. Global quantities, 
however, should still be convergent. A sufficient global 
observable is the gravitational waveform and we perform 
a convergence check on the waveform amplitudes. 

Since we do not have an exact solution to compare 
with, we perform a three- level convergence check, i.e., we 
compute the ratio of the differences in the strain -D/i+,e 
between the three resolutions 

C = Pi ±^ . (74) 

The ratio C defines the convergence factor of the solu- 



tion, and can be translated into the order of convergence 
of the numerics, i.e., the convergence rate. Since our low- 
est order of accuracy is given by first order near shocks, 
we expect at least first-order convergence. 

Checking for convergence in the strain D/i+^e a-s com- 
puted from all extraction methods, we find a convergence 
factor of C > 1 prior to core-bounce and C ~ 1 after 
core bounce. For instance, in Fig. |10[ we show a con- 
vergence plot of the waveform Dh'^^ obtained from the 
CCE method**. In the upper panel, we show the wave- 
forms obtained from three different resolutions, while in 
the lower panel, we show the differences between medium 
and low (blue curve) and high and medium resolutions 
(red curve). Given our resolutions, this convergence fac- 
tor corresponds to a convergence rate between first and 
second order prior to core bounce which reduces to first 
order after core bounce. 

We can estimate a numerical error in the medium res- 
olution simulation by performing a Richardson extrapo- 
lation using the measured convergence rate on the com- 
puted waveform (see e.g. [951). This error estimate, how- 
ever, can only be applied if the convergence rate is un- 
ambiguous. Since we measure a rate between first and 
second order, this is not exactly the case here. Another 
measure of the numerical error in the medium resolution 



* The characteristic computational grid resolutions are scaled by 
the same factors as the corresponding resolutions of the Cauchy 
evolution. 
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FIG. 10: Convergence of the waveform of model A3B3G3 as 
computed at J"^ using CCE. In the bottom panel, we show 
the absolute differences a between low and medium resolu- 
tion, and high and medium resolutions, respectively. With- 
out rescaling any of the difference curves, we observe that they 
approximately line up at and after core bounce so that the 
convergence factor is simply given by C ~ 1. This indicates a 
convergence order at and after core bounce of 1. In the pre- 
bounce phase, the difference between medium and high reso- 
lution is slightly smaller than the difference between medium 
and low resolutions, resulting in a slightly larger convergence 
factor. Given our resolutions, this factor corresponds roughly 
to second-order convergence. 



simulation is therefore directly given by the difference 
between medium and high resolution simulations. This 
is displayed by the green curve in the bottom panel of 
Fig. |10[ At core bounce, we measure in the medium res- 
olution simulation a relative numerical error of ~ 4%. 



V. SUMMARY AND CONCLUSIONS 

We have performed a comparison study of four cur- 
rently available GW extraction techniques in the context 
of axisymmetric rotating stellar core collapse. This study 
is the first to succeed in extracting GWs directly from 
axisymmetric core collapse spacetimes and the first to 
employ the fully coordinate independent CCE extraction 
method for non-vacuum spacetimes. 

We have performed core collapse simulations with sim- 
plified microphysics using a set of three representative 
initial configurations leading to GW signals of varying 
strength and signal morphology in quantitative agree- 
ment with what is expected from microphysically more 
complete models. In addition to having extracted waves 



with variants of the standard coordinate-dependent slow- 
motion, weak-field quadrupole formula, we have em- 
ployed (i) the Regge-Wheeler-Zerilli-Moncrief (RWZM) 
formalism, (ii) extraction based on the Newman-Penrose 
(NP) scalar ^'4, and (iii) Cauchy-characteristic extrac- 
tion (CCE). Of these three latter curvature-based meth- 
ods, RWZM and NP extract GWs at a finite radius from 
the source, and hence, are generally prone to systematic 
errors arising from (i) near-zone effects, (ii) gauge am- 
biguities, and (iii) non-vanishing matter contributions. 
The CCE method, on the other hand, extracts waves 
gauge invariantly at future null infinity J^~^ , that is, at 
an infinite distance from the source where gravitational 
radiation is unambiguously defined. Hence, it is subject 
only to small systematic errors due to the presence of 
matter fields at the CCE world-tube locations. 

An integral ingredient contributing to our success in 
extracting GWs from core collapse using curvature-based 
methods has been the removal of unphysical non-linear 
low-frequency drifts from the waveforms that otherwise 
would make a proper analysis largely impossible. This 
has been achieved by the application of fixed-frequency 
integration (FFI, [102j ) for time integration and filtering 
to yield the strain h. 

Comparing the waveforms obtained with the various 
extraction methods, we make a number of observations: 

(i) NP- and CCE-extracted waveforms converge with ex- 
traction and world-tube radius, respectively. The wave- 
forms obtained with the RWZM formalism show spurious 
high-frequency components that no other method repro- 
duces. A number of tests imply that the RWZM method 
may be less applicable to weak GW signals, at least at the 
currently accessible numerical resolutions and grid sizes. 

(ii) NP extraction, CCE, and even the quadrupole ap- 
proximation, yield waveforms which agree well in phase, 
with differences in the time lags between successive peaks 
of < 0.05 ms. Since the RWZM formalism is contami- 
nated by unphysical high-frequency components, an ac- 
curate determination of the phasing compared to the 
other methods is largely impossible, (iii) The maximum 
amplitudes at core bounce are different by ~ 1 — 7% in 
waveforms obtained with NP extraction and are system- 
atically smaller by '-^ 5 — 11% in waveforms obtained with 
the QF compared to the waves obtained via CCE. Ac- 
cordingly, CCE yields waveforms that result in slightly 
higher signal-to-noise-ratios (SNRs) (^ 6 — 9%). (iv) 
Overall, the error of the waveforms computed with the 
quadrupole approximation are well within numerical er- 
rors and physical uncertainties. Unlike the waveforms ob- 
tained with the curvature-based methods, the quadrupole 
waveforms do not suffer from low-frequency drifts. In 
that respect, the quadrupole approximation is advanta- 
geous. We also observe that the quadrupole variant using 
"physical" velocity components [75] yields waves that are 
closer to those obtained via CCE. However, this finding 
may be true only for the core collapse case studied here 
and may not hold in general, (v) While it is unlikely that 
matched filtering approaches will be used in searches for 
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GWs from core collapse in the near future, we have never- 
theless computed GW template mismatches, a measure 
for the detectability of differences between waveforms. 
We find that when used in hypothetical matched-filtering 
GW searches, waveforms from NP extraction, CCE, and 
the QF would lead to the detection of the same model, 
while the waveforms computed with the RWZM formal- 
ism would generally not. 

There are two major drawbacks of our current work: 
(i) The curvature-based methods assume vacuum at the 
extraction spheres and world-tube locations. Hence, we 
must, in principle, extract at very large radii where the 
stress-energy tensor is zero. This, however, is currently 
not possible, since the collapsing star extends over the 
entire computational grid and larger grids are compu- 
tationally prohibitive, (ii) All curvature-based methods 
yield waveforms with unphysical low-frequency drifts, re- 
quiring removal by spectral cut-off via FFI. This is partic- 
ularly problematic in models with physical content below 
~ 100 Hz. A possible improvement of the low-frequency 
behavior could be achieved by the inclusion of matter 
terms in the CCE method, or alternatively, by enlarg- 
ing the simulation domain such that the extraction takes 
place outside of the star and in pure vacuum. The lat- 
ter could be efficiently achieved by employing multiblock 
techniques that cover the wavezone by a set of spherical 
grids jll6| . 

Finally, we point out that we have considered only the 
GW signal from rotating core collapse and bounce in this 
first study using curvature-based GW extraction from 
core collapse spacetimes. While our results may transfer 
to other GW emission processes in core collapse, this is 
by no means guaranteed. Further work will be needed 
to adress curvature-based GW extraction also from post- 
bounce convection and the standing accretion shock in- 



stability, protoneutron star pulsations, rotational insta- 
bilities, and black hole formation. 
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